11 Image data

## ----data1, fig.cap=""--------------------------------------------------------
library("EBImage")
imagefile = system.file("images", "mosquito.png",
                        package = "MSMB")
mosq = readImage(imagefile)
## ----vis1, eval = FALSE-------------------------------------------------------
display(mosq)

## ----mosquito, fig.keep = 'high', fig.cap = "Mosquito discovered deceased in the suburbs of Decatur, Georgia (credit: CDC / Janice Haney Carr).", dev = "png", dpi = 300, fig.width=dim(mosq)[1]/300, fig.height=dim(mosq)[2]/300----
display(mosq, method = "raster")
text(x = 85, y = 800, label = "A mosquito",
     adj = 0, col = "orange", cex = 1.5)

display(1-mosq, method = "raster")

## ----vis3a, eval = FALSE------------------------------------------------------
imagefile = system.file("images", "hiv.png",
                        package = "MSMB")
hivc = readImage(imagefile)
display(hivc)

## ----vis3b, eval = TRUE, echo = FALSE-----------------------------------------
imagefile = system.file("images", "hiv.png",
                        package = "MSMB")
hivc = readImage(imagefile)

## ---- hiv, eval = TRUE, echo = FALSE, fig.show = 'hold', fig.keep = 'high', fig.cap = "Scanning electron micrograph of HIV-1 virions budding from a cultured lymphocyte (credit: CDC / C. Goldsmith, P. Feorino, E.L. Palmer, W.R. McManus)."----
knitr::include_graphics(c('../images/hiv.png'), dpi = NA)

## ----image-oneminus, fig.keep = 'high', fig.cap = "Tiled display of four images of cell nuclei from the **[EBImage](https://bioconductor.org/packages/EBImage/)** package.", fig.margin = FALSE, dev = "png"----
nuc = readImage(system.file("images", "nuclei.tif",
                            package = "EBImage"))

nuc
## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 510 510 4 
##   frames.total : 4 
##   frames.render: 4 
## 
## imageData(object)[1:5,1:6,1]
##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,] 0.06274510 0.07450980 0.07058824 0.08235294 0.10588235 0.09803922
## [2,] 0.06274510 0.05882353 0.07843137 0.09019608 0.09019608 0.10588235
## [3,] 0.06666667 0.06666667 0.08235294 0.07843137 0.09411765 0.09411765
## [4,] 0.06666667 0.06666667 0.07058824 0.08627451 0.08627451 0.09803922
## [5,] 0.05882353 0.06666667 0.07058824 0.08235294 0.09411765 0.10588235
dim(imageData(nuc))
## [1] 510 510   4
display(1 - nuc, method = "raster", all = TRUE)

display( nuc, method = "raster", all = TRUE)

## ----nucleitiledoneframe, dev = "png"-----------------------------------------
display(1 - nuc, method = "raster", frame = 2)

## ----how1---------------------------------------------------------------------
class(mosq)
## [1] "Image"
## attr(,"package")
## [1] "EBImage"
mosq[1:3,1:3]
## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 3 3 
##   frames.total : 1 
##   frames.render: 1 
## 
## imageData(object)[1:3,1:3]
##           [,1]      [,2]      [,3]
## [1,] 0.1960784 0.1960784 0.1960784
## [2,] 0.1960784 0.1960784 0.1960784
## [3,] 0.1960784 0.1960784 0.2000000
## ----how2, echo = FALSE, eval = FALSE-----------------------------------------
showClass("Image")
## Class "Image" [package "EBImage"]
## 
## Slots:
##                           
## Name:      .Data colormode
## Class:     array   integer
## 
## Extends: 
## Class "array", from data part
## Class "structure", by class "array", distance 2
## Class "vector", by class "array", distance 3, with explicit coerce
## ----dim----------------------------------------------------------------------
dim(mosq)
## [1] 1400  952
## ----mosqhist, fig.keep = 'high', fig.cap = "Histogram of the pixel intensities in `mosq`. Note that the range is between 0 and 1.", fig.width = .h$width, fig.height = .h$height----
hist(mosq)

## ----check, echo=FALSE--------------------------------------------------------
stopifnot(all(mosq>=0 & mosq<=1), isTRUE(all.equal(max(mosq), 1)), isTRUE(all.equal(min(mosq), 0)))


## ----how3---------------------------------------------------------------------
dim(imageData(mosq))
## [1] 1400  952
imageData(mosq)[1:3, 1:6]
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784
## [2,] 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784
## [3,] 0.1960784 0.1960784 0.2000000 0.2039216 0.2000000 0.1960784
## ----show1--------------------------------------------------------------------
mosq
## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 1400 952 
##   frames.total : 1 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6]
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784
## [2,] 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784
## [3,] 0.1960784 0.1960784 0.2000000 0.2039216 0.2000000 0.1960784
## [4,] 0.1960784 0.1960784 0.2039216 0.2078431 0.2000000 0.1960784
## [5,] 0.1960784 0.2000000 0.2117647 0.2156863 0.2000000 0.1921569
## ----show2--------------------------------------------------------------------
hivc
## Image 
##   colorMode    : Color 
##   storage.mode : double 
##   dim          : 1400 930 3 
##   frames.total : 3 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6,1]
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    0    0    0    0    0    0
## [2,]    0    0    0    0    0    0
## [3,]    0    0    0    0    0    0
## [4,]    0    0    0    0    0    0
## [5,]    0    0    0    0    0    0
## ----show3, echo=FALSE--------------------------------------------------------
stopifnot(colorMode(mosq)==Grayscale, colorMode(hivc)==Color, dim(nuc)[3]==4)


## ----how6, results = "hide"---------------------------------------------------
nuc
## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 510 510 4 
##   frames.total : 4 
##   frames.render: 4 
## 
## imageData(object)[1:5,1:6,1]
##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,] 0.06274510 0.07450980 0.07058824 0.08235294 0.10588235 0.09803922
## [2,] 0.06274510 0.05882353 0.07843137 0.09019608 0.09019608 0.10588235
## [3,] 0.06666667 0.06666667 0.08235294 0.07843137 0.09411765 0.09411765
## [4,] 0.06666667 0.06666667 0.07058824 0.08627451 0.08627451 0.09803922
## [5,] 0.05882353 0.06666667 0.07058824 0.08235294 0.09411765 0.10588235
dim(imageData(nuc))
## [1] 510 510   4
## ----checkassertion, echo = FALSE---------------------------------------------
stopifnot(all(c("  frames.total : 4 ", "  frames.render: 4 ") %in%
              capture.output(EBImage:::showImage(nuc))))


## ----write1-------------------------------------------------------------------
writeImage(hivc, "hivc.jpeg", quality = 85)
dim(hivc)
## [1] 1400  930    3
## ----objectsize---------------------------------------------------------------
object.size(hivc) |> format(units = "Mb")
## [1] "29.8 Mb"
object.size(hivc) |> format(units = "Kb")
## [1] "30516.4 Kb"
(object.size(hivc) / prod(dim(hivc))) |> format() |> paste("per pixel")
## [1] "8 bytes per pixel"
file.info("hivc.jpeg")$size
## [1] 294904
16 * 3 * 8 #three color, 16 Megapixel image
## [1] 384
## ---- mosqcrop, eval = TRUE, echo = FALSE, fig.show = 'hold', fig.keep = 'high', fig.cap = "Cropping: `mosqcrop`"----
knitr::include_graphics(c('../images/mosquito.png'), dpi = NA)

## ----manip1a------------------------------------------------------------------
mosqinv = normalize(-mosq)


## ----manip3a------------------------------------------------------------------
mosqcont = mosq * 3
mosqexp = mosq ^ (1/3)
display(mosqcont, method = "raster")

display(mosqexp, method = "raster")

display(mosq+0.4, method = "raster")

display(transpose(mosq), method = "raster")

mosq[100:438, 112:550]
## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 339 439 
##   frames.total : 1 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6]
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.1882353 0.1882353 0.1921569 0.1882353 0.1882353 0.1882353
## [2,] 0.1921569 0.1921569 0.1921569 0.1960784 0.2000000 0.1921569
## [3,] 0.1960784 0.1960784 0.1960784 0.2078431 0.2078431 0.1960784
## [4,] 0.1960784 0.1882353 0.1843137 0.1960784 0.2039216 0.2039216
## [5,] 0.1960784 0.1803922 0.1764706 0.1882353 0.2000000 0.2078431
## ----manip4a------------------------------------------------------------------
mosqcrop   = mosq[100:438, 112:550]
mosqthresh = mosq > 0.5
mosqtransp = transpose(mosq)
display(mosqcrop, method = "raster")

display(mosqthresh, method = "raster")

display(mosqtransp, method = "raster")

## ----checkassertionont, echo = FALSE------------------------------------------
stopifnot(identical(t(mosq), transpose(mosq)))
## ----spattrans1---------------------------------------------------------------
mosqrot   = EBImage::rotate(mosq, angle = 30)
mosqshift = translate(mosq, v = c(40, 70))
mosqflip  = flip(mosq)
mosqflop  = flop(mosq)

display(mosqrot, method = "raster")

display(mosqshift, method = "raster")

display(mosqflip, method = "raster")

display(mosqflop, method = "raster")

## ----MSMB, results="hide"-----------------------------------------------------
imagefiles = system.file("images", c("image-DAPI.tif",
  "image-FITC.tif", "image-Cy3.tif"), package="MSMB")
cells = readImage(imagefiles)
display(cells)
## Only the first frame of the image stack is displayed.
## To display all frames use 'all = TRUE'.

## ----checkdim, echo=FALSE-----------------------------------------------------
stopifnot(dim(cells)[3]==3)

## ----range--------------------------------------------------------------------
apply(cells, 3, range)
##       image-DAPI  image-FITC   image-Cy3
## [1,] 0.001586938 0.002899214 0.001663233
## [2,] 0.031204700 0.062485695 0.055710689
## ----fixrange-----------------------------------------------------------------
cells[,,1]   = 32 * cells[,,1]
cells[,,2:3] = 16 * cells[,,2:3]
apply(cells, 3, range)
##      image-DAPI image-FITC  image-Cy3
## [1,] 0.05078202 0.04638743 0.02661173
## [2,] 0.99855039 0.99977111 0.89137102
getFrame(cells, 1)
## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 340 490 
##   frames.total : 1 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6]
##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,] 0.05517662 0.05420005 0.05566491 0.05371176 0.05517662 0.05322347
## [2,] 0.05566491 0.05517662 0.05420005 0.05420005 0.05566491 0.05468833
## [3,] 0.05371176 0.05615320 0.05517662 0.05517662 0.05371176 0.05322347
## [4,] 0.05517662 0.05468833 0.05420005 0.05566491 0.05420005 0.05273518
## [5,] 0.05566491 0.05468833 0.05664149 0.05371176 0.05517662 0.05371176
## ----defw---------------------------------------------------------------------
w = makeBrush(size = 51, shape = "gaussian", sigma = 7)
nucSmooth = filter2(getFrame(cells, 1), w)
dim(w)
## [1] 51 51
## ----image-filter2, fig.keep = 'high', fig.cap = "The middle row of the weight matrix, `w[` (ref:image-filter2-1) `, ]`", fig.width = 3, fig.height = 2.6, dev = "png"----
library("tibble")
library("ggplot2")
tibble(w = w[(nrow(w)+1)/2, ]) |>
  ggplot(aes(y = w, x = seq(along = w))) + geom_point()

## ----smooth-------------------------------------------------------------------
cellsSmooth = Image(dim = dim(cells))
sigma = c(1, 3, 3)
for(i in seq_along(sigma))
  cellsSmooth[,,i] = filter2( cells[,,i],
         filter = makeBrush(size = 51, shape = "gaussian",
                            sigma = sigma[i]) )
## ----illuminationartifact1----------------------------------------------------
py = seq(-1, +1, length.out = dim(cellsSmooth)[1])
px = seq(-1, +1, length.out = dim(cellsSmooth)[2])
illuminationGradient = Image(
     outer(py, px, function(x, y) exp(-(x^2+y^2))))
nucBadlyIlluminated = cellsSmooth[,,1] * illuminationGradient
## ----illuminationartifact2----------------------------------------------------
disc = makeBrush(21, "disc")
disc = disc / sum(disc)
localBackground = filter2(nucBadlyIlluminated, disc)
offset = 0.02
nucBadThresh = (nucBadlyIlluminated - localBackground > offset)
## ----adathresh----------------------------------------------------------------
nucThresh =
  (cellsSmooth[,,1] - filter2(cellsSmooth[,,1], disc) > offset)


## ----morphopen1---------------------------------------------------------------
nucOpened = EBImage::opening(nucThresh,
                  kern = makeBrush(5, shape = "disc"))
## ----imageProcessing14--------------------------------------------------------
nucSeed = bwlabel(nucOpened)
table(nucSeed)
## nucSeed
##      0      1      2      3      4      5      6      7      8      9     10 
## 155408    511    330    120    468    222    121    125    159    116    520 
##     11     12     13     14     15     16     17     18     19     20     21 
##    115    184    179    116    183    187    303    226    164    309    194 
##     22     23     24     25     26     27     28     29     30     31     32 
##    148    345    287    203    379    371    208    222    320    443    409 
##     33     34     35     36     37     38     39     40     41     42     43 
##    493    256    169    225    376    214    228    341    269    119    315
## ----imageProcessing17, eval = FALSE------------------------------------------
## display(colorLabels(nucSeed))

## ----imageProcessing15a-------------------------------------------------------
nucMask = cellsSmooth[,,1] - filter2(cellsSmooth[,,1], disc) > 0


## ----imageProcessing15b-------------------------------------------------------
nucMask = fillHull(nucMask)


## ----imageProcessing16--------------------------------------------------------
nuclei = propagate(cellsSmooth[,,1], nucSeed, mask = nucMask)
## ----voronoiExample-----------------------------------------------------------
zeros        = Image(dim = dim(nuclei))
voronoiExamp = propagate(seeds = nuclei, x = zeros, lambda = 100)
voronoiPaint = paintObjects(voronoiExamp, 1 - nucOpened)


## ----voronoiEx----------------------------------------------------------------
head(table(voronoiExamp))
## voronoiExamp
##    1    2    3    4    5    6 
## 5645 4735  370 5964 3333 1377
ind = which(voronoiExamp == 13, arr.ind = TRUE)
head(ind, 3)
##      row col
## [1,] 112 100
## [2,] 113 100
## [3,] 114 100
## ----histcellbody-1, fig.width = .h$width, fig.height = .h$height, fig.keep="high", eval = TRUE, echo = FALSE, fig.cap = "Histogram of the actin channel in `cellsSmooth`, after taking the logarithm."----
hist(log(cellsSmooth[,,3]) )

## ----histcellbody-2, fig.width = .h$width, fig.height = .h$height, fig.keep="high", eval = TRUE, echo = FALSE, fig.cap = "Zoom into Figure \\@ref(fig:histcellbody-1)."----
hist(log(cellsSmooth[,,3]), xlim = -c(3.6, 3.1), breaks = 300)

## ----histcellbody, eval = FALSE, echo = TRUE----------------------------------
## hist(log(cellsSmooth[,,3]) )
## hist(log(cellsSmooth[,,3]), xlim = -c(3.6, 3.1), breaks = 300)


## ----checkhistcellsSmooth, echo=FALSE-----------------------------------------
stopifnot(mean(cellsSmooth[,,3]>=exp(-3.6) & cellsSmooth[,,3]<=exp(-3.1)) > 0.68)


## ----musigmaEstimator, message=FALSE------------------------------------------
library("genefilter")
bgPars = function(x) {
  x    = log(x)
  loc  = half.range.mode( x )
  left = (x - loc)[ x < loc ]
  wid  = sqrt( mean(left^2) )
  c(loc = loc, wid = wid, thr = loc + 6*wid)
}
cellBg = apply(cellsSmooth, MARGIN = 3, FUN = bgPars)
cellBg
##            [,1]        [,2]        [,3]
## loc -2.90176965 -2.94427499 -3.52191681
## wid  0.00635322  0.01121337  0.01528207
## thr -2.86365033 -2.87699477 -3.43022437
## ----histcellbody-3, fig.keep = 'high', fig.cap = "As in Figure \\@ref(fig:histcellbody-2), but with `loc` and `thr` shown by vertical lines.", fig.width = .h$width, fig.height = .h$height----
hist(log(cellsSmooth[,,3]), xlim = -c(3.6, 3.1), breaks = 300)
abline(v = cellBg[c("loc", "thr"), 3], col = c("brown", "red"))

## ----cytoplasmMask------------------------------------------------------------
cytoplasmMask = (cellsSmooth[,,2] > exp(cellBg["thr", 2])) |
       nuclei | (cellsSmooth[,,3] > exp(cellBg["thr", 3]))


## ----imageProcessing22--------------------------------------------------------
cellbodies = propagate(x = cellsSmooth[,,3], seeds = nuclei,
                       lambda = 1.0e-2, mask = cytoplasmMask)
## ----imageProcessing25--------------------------------------------------------
cellsColor = rgbImage(red   = cells[,,3],
                      green = cells[,,2],
                      blue  = cells[,,1])

nucSegOnNuc  = paintObjects(nuclei, tgt = toRGB(cells[,,1]),
                            col = "#ffff00")
nucSegOnAll  = paintObjects(nuclei, tgt = cellsColor,
                            col = "#ffff00")
cellSegOnAll = paintObjects(cellbodies, tgt = nucSegOnAll,
                            col = "#ff0080")
## ----baserfeats---------------------------------------------------------------
meanNucInt       = tapply(cells[,,1], nuclei, mean)
meanActIntInNuc  = tapply(cells[,,3], nuclei, mean)
meanActIntInCell = tapply(cells[,,3], cellbodies, mean)


## ----pairsint, fig.keep = 'high', fig.cap = "Pairwise scatterplots of per-cell intensity descriptors. \\label{pairsint}", fig.margin = FALSE, fig.width = 5.3, fig.height = 5.3----
library("GGally")
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(tibble(meanNucInt, meanActIntInNuc, meanActIntInCell))

## ----imageProcessing27--------------------------------------------------------
F1 = computeFeatures(nuclei,     cells[,,1], xname = "nuc",
                                             refnames = "nuc")
F2 = computeFeatures(cellbodies, cells[,,2], xname = "cell",
                                             refnames = "tub")
F3 = computeFeatures(cellbodies, cells[,,3], xname = "cell",
                                             refnames = "act")
dim(F1)
## [1] 43 89
## ----showF1-------------------------------------------------------------------
 F1[1:3, 1:5]
##   nuc.0.m.cx nuc.0.m.cy nuc.0.m.majoraxis nuc.0.m.eccentricity nuc.0.m.theta
## 1   119.5523   17.46895          44.86819            0.8372059     -1.314789
## 2   143.4511   15.83709          26.15009            0.6627672     -1.213444
## 3   336.5401   11.48175          18.97424            0.8564444      1.470913
## ----readlymphnodedata--------------------------------------------------------
library("readr")
## 
## Attaching package: 'readr'
## The following object is masked from 'package:genefilter':
## 
##     spec
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:EBImage':
## 
##     combine
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
cellclasses = c("T_cells", "Tumor", "DCs", "other_cells")
brcalymphnode = lapply(cellclasses, function(k) {
    read_csv(file.path("..", "data",
             sprintf("99_4525D-%s.txt", k))) |>
    transmute(x = globalX,
              y = globalY,
              class = k)
}) |> bind_rows() |> mutate(class = factor(class))
## Rows: 103681 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): filename
## dbl (4): locX, locY, globalX, globalY
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 27822 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): filename
## dbl (4): locX, locY, globalX, globalY
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 878 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): filename
## dbl (4): locX, locY, globalX, globalY
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 77081 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): filename
## dbl (4): locX, locY, globalX, globalY
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
brcalymphnode
## # A tibble: 209,462 × 3
##        x     y class  
##    <dbl> <dbl> <fct>  
##  1  6355 10382 T_cells
##  2  6356 10850 T_cells
##  3  6357 11070 T_cells
##  4  6357 11082 T_cells
##  5  6358 10600 T_cells
##  6  6361 10301 T_cells
##  7  6369 10309 T_cells
##  8  6374 10395 T_cells
##  9  6377 10448 T_cells
## 10  6379 10279 T_cells
## # ℹ 209,452 more rows
table(brcalymphnode$class)
## 
##         DCs other_cells     T_cells       Tumor 
##         878       77081      103681       27822
## ----checkcellnumber, echo = FALSE--------------------------------------------
tabtab = table(brcalymphnode$class)
within = function(x, a, b) (x>a & x<b)
stopifnot(all(within(tabtab[c("T_cells", "Tumor", "DCs")], c(100000, 27000, 800), c(110000, 28000, 1000))))


## ----brcalntcells, fig.keep = 'high', fig.cap = "Scatterplot of the $x$ and $y$ positions of the T- and tumor cells in `brcalymphnode`. The locations were obtained by a segmentation algorithm from a high resolution version of Figure \\@ref(fig:stainedlymphnode). Some rectangular areas in the T-cells plot are suspiciously empty, this could be because the corresponding image tiles within the overall composite image went missing, or were not analyzed.", fig.margin = FALSE, dev = "png", dpi = 300, fig.width=9, fig.height=4.5, pointsize=24----
ggplot(filter(brcalymphnode, class %in% c("T_cells", "Tumor")),
   aes(x = x, y = y, col = class)) + geom_point(shape = ".") +
   facet_grid( . ~ class) + guides(col = "none")

## ----spatstat1, results = "hide"----------------------------------------------
library("spatstat")
## Loading required package: spatstat.data
## 
## Attaching package: 'spatstat.data'
## The following object is masked _by_ '.GlobalEnv':
## 
##     cells
## Loading required package: spatstat.geom
## spatstat.geom 3.0-3
## 
## Attaching package: 'spatstat.geom'
## The following object is masked from 'package:genefilter':
## 
##     area
## The following objects are masked from 'package:EBImage':
## 
##     affine, closing, distmap, opening, rotate
## Loading required package: spatstat.random
## spatstat.random 3.0-1
## Loading required package: spatstat.core
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Loading required package: rpart
## spatstat.core 2.4-0
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-2
## 
## spatstat 2.3-3       (nickname: 'That's not important right now') 
## For an introduction to spatstat, type 'beginner'
## ----spatstat2----------------------------------------------------------------
ln = with(brcalymphnode,
  ppp(x = x, y = y, marks = class, xrange = range(x), yrange = range(y)))
ln
## Marked planar point pattern: 209462 points
## Multitype, with levels = DCs, other_cells, T_cells, Tumor 
## window: rectangle = [3839, 17276] x [6713, 23006] units
## ----checkclassln, echo = FALSE-----------------------------------------------
stopifnot(identical(class(ln), "ppp"))
## ----spatstat3----------------------------------------------------------------
library("spatstat")
ln = with(brcalymphnode, ppp(x = x, y = y, marks = class, 
                             xrange = range(x), yrange = range(y)))
ln
## Marked planar point pattern: 209462 points
## Multitype, with levels = DCs, other_cells, T_cells, Tumor 
## window: rectangle = [3839, 17276] x [6713, 23006] units
## ----densityppz1--------------------------------------------------------------
d = density(subset(ln, marks == "Tumor"), edge=TRUE, diggle=TRUE)
plot(d)

## ----densityppp1, fig.keep = 'high', fig.cap = "Intensity estimate for the cells marked `Tumor` in `ppp`. The support of the estimate is the polygon that we specified earlier on (Figure \\@ref(fig:convhull)).", fig.width = 3.75, fig.height = 3.5, echo = FALSE----
par(mai = c(0, 0, 0.2, 0.7))
plot(d)

## ----densityppz0--------------------------------------------------------------
d0 = density(subset(ln, marks == "Tumor"), edge = FALSE)
plot(d0)

## ----densityppp0, fig.keep = 'high', fig.cap = "As Figure \\@ref(fig:densityppp1), but without edge correction \\label{densityppp0}", fig.width = 3.75, fig.height = 3.5, echo = FALSE----
par(mai = c(0, 0, 0.2, 0.7))
plot(d0)

## ----relrisk-calc-------------------------------------------------------------
rr = relrisk(ln, sigma = 250)


## ----relrisk, fig.keep = 'high', fig.cap = "Estimates of the spatially varying probability of each of the cell claases, conditional on there being cells.", fig.margin = FALSE, fig.width=7, fig.height=7----
plot(rr)

## ----checkConditional, echo=FALSE---------------------------------------------
m = rr[[1]]$v
for(i in 2:length(rr)) m = m + rr[[i]]$v
#stopifnot(all(is.na(m) | abs(m-1)<1e-6))


## ----Gestshow-----------------------------------------------------------------
gln = Gest(ln)
gln
## Function value object (class 'fv')
## for the function r -> G(r)
## .....................................................................
##         Math.label      Description                                  
## r       r               distance argument r                          
## theo    G[pois](r)      theoretical Poisson G(r)                     
## han     hat(G)[han](r)  Hanisch estimate of G(r)                     
## rs      hat(G)[bord](r) border corrected estimate of G(r)            
## km      hat(G)[km](r)   Kaplan-Meier estimate of G(r)                
## hazard  hat(h)[km](r)   Kaplan-Meier estimate of hazard function h(r)
## theohaz h[pois](r)      theoretical Poisson hazard function h(r)     
## .....................................................................
## Default plot formula:  .~r
## where "." stands for 'km', 'rs', 'han', 'theo'
## Recommended range of argument r: [0, 21.033]
## Available range of argument r: [0, 61.89]
## ----Gest, fig.keep = 'high', fig.cap = "Estimates of $G$, using three different edge effect corrections --which here happen to essentially lie on top of each other-- and the theoretical distribution for a homogenous Poisson process. \\label{Gest}", fig.width = 4.25, fig.height = 4.25, results="hide"----
library("RColorBrewer")
plot(gln, xlim = c(0, 10), lty = 1, col = brewer.pal(4, "Set1"))

## ----Linhom-------------------------------------------------------------------
Lln = Linhom(subset(ln, marks == "T_cells"))
## number of data points exceeds 1000 - computing border correction estimate only
Lln
## Function value object (class 'fv')
## for the function r -> L[inhom](r)
## ................................................................................
##            Math.label                
## r          r                         
## theo       L[pois](r)                
## border     {hat(L)[inhom]^{bord}}(r) 
## bord.modif {hat(L)[inhom]^{bordm}}(r)
##            Description                                      
## r          distance argument r                              
## theo       theoretical Poisson L[inhom](r)                  
## border     border-corrected estimate of L[inhom](r)         
## bord.modif modified border-corrected estimate of L[inhom](r)
## ................................................................................
## Default plot formula:  .~.x
## where "." stands for 'bord.modif', 'border', 'theo'
## Recommended range of argument r: [0, 819.84]
## Available range of argument r: [0, 819.84]
## ----Images-Lln, fig.keep = 'high', fig.cap = "Estimate of $L_{\\scriptsize \\mbox{inhom}}$, Equations \\@ref(eq:kinhom) and \\@ref(eq:Lest), of the T cell pattern. \\label{Images-Lln}", fig.width = 4, fig.height = 5.9, results = "hide"----
plot(Lln, lty = 1, col = brewer.pal(3, "Set1"))

## ----pcfdo--------------------------------------------------------------------
pcfln = pcf(Kinhom(subset(ln, marks == "T_cells")))
## number of data points exceeds 1000 - computing border correction estimate only
## ----Images-pcf, fig.keep = 'high', fig.show = "hold", fig.cap = "Estimate of the pair correlation function, Equation \\@ref(eq:pcf), of the T cell pattern. \\label{Images-pcf}", fig.width=5, fig.height=5, results="hide"----
plot(pcfln, lty = 1)

plot(pcfln, lty = 1, xlim = c(0, 10))

## ----samplingpcf, fig.keep = 'high', fig.cap = "Answer to Question \\@ref(ques:Image-ques-samplingpcf): as in the bottom panel of Figure \\@ref(fig:Images-pcf), but with denser sampling.", fig.width=5, fig.height=5----
pcfln2 = pcf(Kinhom(subset(ln, marks == "T_cells"), r = seq(0, 10, by = 0.2)))
## number of data points exceeds 1000 - computing border correction estimate only
plot(pcfln2, lty = 1)